  function [G1,C,HH,eu,SDX,zmat,NY,NX,NETA,gev,aux]=play_eugen(modstr,param); 
% function [G1,C,HH,eu,SDX,zmat,NY,NX,NETA,gev,aux]=play_eugen(modstr,param); 
%
% This function makes a call to gensys for the model in modstr witha a given 
% set of coefficients called in param and spits out the determinacy 
% eigenvalue count, and checks dimensions of the matrices 
% Error will occur if there are dimension mistmaches or indeterminacy 
% 
% Alejandro Justiniano Nov 24 2006 
% =======================================================================
[G1,C,HH,eu,SDX,zmat,NY,NX,NETA,gev]=feval(modstr,param);

disp(' '); 

% =======================================
aux.ny=NY;aux.nx=NX;aux.nz=size(zmat,1); 


if isstr( modstr ) == 0 
    modstr = func2str( modstr ); 
end 

disp(['Begin Analysis for model: ',modstr]); 

%Compare lengths 
%===============
ng=size(G1);
if ~isempty(G1)
    if ng(1)~=NY | ng(2) ~= NY;
        error('Mistmatch NY or G1')
    end
else 
    disp('G1 is empty, model allows indeterminacy?'); 
    quer('c') 
    return
end 

if ~isempty(HH);
    nimp=size(HH);
    if nimp(1) ~= NY | nimp(2)~=NX
        error('Mistmatch NX or HH');
    end 
else 
    disp('HH is empty, model allows indeterminacy?');
    quer('c')
    return
end 
%if any( size(SDX) ~= [NX NX] == 1); 
%    error('Mistmatch SDX and NX'); 
%end; 

disp(' '); 
disp('Model Dimensions'); 
dispaj('Number of NY variables:',NY); 
dispaj('Number of Shocks variables:',NX); 

% =================
% Check Zmat matrix 
% =================
nz=size(zmat);
if nz(2) ~= ng(1);
    error('Size of zmat matrix does not match G1')
end;

disp(' '); 
dispaj('Number of Observable Series :',nz(1) ); 
temp=find( all( zmat' == 0 ) == 1 ); 
if ~isempty( temp ); 
    disp('All zeros in following rows of zmat matrix');
    disp( num2str( temp ) ); 
    error('All zeros in rows of zmat matrix'); 
end 

temp=zmat( find(zmat~=0) ); 
disp('Coefficients in zmat'); 
disp( num2str(temp) ); 

% Analyze gev 
if ~isempty(gev)
    % ===========
    gev2 = gev(:,2) ./ gev(:,1);
    ns=length( find(gev2 < 0.9999 ) );
    nu=length( find(gev2 > 1.0001 ) );
    ng=length(gev2)  ;
    nuint=length( find( abs(gev2 - ones(ng,1)) < eps*1000 ) );

    if isempty(nuint)==1 ;
        nuint = 0 ;
    end
    disp('Eigenvalue count');
    eigcell=cell(3,3);
    eigcell(1,:)={'Stable','Unstable','Unit'};
    eigcell(2,:)={'======','========','====='};
    eigcell(3,:)=num2cellaj([ns nu nuint]);
    printcell(eigcell);

    if eu(1) == 0 | eu(1) == -2
        error( ' Solution does not exist ');
    else
        if eu(2) == 1
            disp('Unique solution exists');
        elseif eu(2) == 0;
            error(' Indeterminacy ')
        end
    end
else 
    disp('GEV empty, solution muste be indeterminate') 
end 

